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Abstract 

The explicit quasi-monotonic conservative TVD scheme and numerical method for 
the solution of the gravitational MHD equations are developed. The 2D numerical 
code for the simulation of multidimensional selfgravitating MHD flows on the Eu- 
lerian cylindrical grid is constructed. The results of test calculations show that the 
code has a good mathematical and computational properties and can be applied to 
the solution of a wide class of plasma physics and astrophysics problems. We have 
simulated, for example, a collapse of magnetized rotating protostellar clouds. 
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1 Introduction 



Magnetic field plays an important role in the various problems of plasma physics and 
astrophysics. It is well known, for example, that the density of the magnetic energy in the 
interstellar clouds is compared with the density of the turbulent energy and can exceed the 
density of the thermal energy (see, e.g., [1]). There are a wide class of problems dealing with 
the dynamics of the selfgravitating MHD flows. Note among others the gravitational collapse 
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of the interstellar clouds, protostellar clusters and galaxies; origin of the molecular bipolar 
outflows and jets from young stars and galaxy cores; dynamics of the accretion flows and 
circumstellar disks. 

A variety of finite-difference schemes for the solution of the gasdynamical equations have 
been developed during last five decades. Corresponding ID, 2D, and 3D numerical codes have 
been applied to the simulation of various gasdynamical flows. The direct extension of the 
these schemes to the MHD, however, is quite complex problem due to the specific properties 
of the MHD equations (anisotropy, vanishing divergence of the magnetic field, etc.). 

Among the non-monotonic generalization of gasdynamical schemes for numerical simulation 
of the MHD fiows the Lax-Vendroff method (and its modifications) is used (see, e.g., [2]). This 
method is subjected to the non-physical oscillations on strong discontinuities [3]. Another 
prevailing method for the simulation of the multidimensional selfgravitating MHD fiows is 
the MHD extension of the SPH method [4]. 

A finite-difference scheme should approximate correctly the physical oscillations of the 
fiows containing shock waves, contact discontinuities and rarefaction waves. The nonlinear 
monotonic methods of high resolution were suggested for the solution of these problems. 
Historically the first monotonic method was FCT method [5,6]. It was developed for of 
linear advection equation and for gasdynamical equations. At the present time FCT schemes 
are used for the simulations of the MHD fiows as well [7] . 

There are two basic approaches that are applied for the construction of the monotonic 
finite-difference schemes for the numerical solution of the MHD problems. Nonlinear Riemann 
problem is solved in the original Godunov's methods at every mesh interface and every time 
step by various iteration procedures (see, e.g., [8]). On the basis of nonlinear Riemann solver, 
Dai and Woodward [9] extended the PPM method [10] there. Powell et al. [11,12] constructed 
an eight-wave eigensystem for the approximate Riemann solver. Jiang and Wu [13] applied 
a high-order WENO interpolation scheme to MHD equations. Another approach is splitting 
the hyperbolicity matrix onto positive and negative parts and following solution of the linear 
Riemann problem [14-20]. 

We have suggested recently a simple TVD scheme for the numerical solution of the MHD 
equations that has high resolution in the smooth regions of the solution (see [21]). This 
scheme does not demand the solution of full eigenvalues problem (i.e. a calculation of all 
eigenvectors and eigenvalues of the hyperbolicity matrix), since it uses only the spectral 
radius. On the basis of developed TVD scheme we elaborated 2D numerical method for the 
solution of the selfgravitational MHD equations in cylindrical coordinates [22]. This method 
can be used for the simulation of isothermal and adiabatic plasma fiows. We have applied it 
to simulation of the collapse of magnetized rotating protostellar clouds. 

In outline this paper proceeds as follows. In Section 2 we consider the mathematical problem 
statement for the simulation of selfgravitating MHD fiows in the frame of axisymmetrical 
approximation. In Section 3 the TVD scheme for the numerical solution of the ID MHD 
system is constructed. In Section 4 the developed ID TVD scheme is extended to the 2D 
case and numerical algorithm for simulation of the 2D axisymmetrical selfgravitating MHD 
fiows is described. Section 5 gives the description of 2D MHD numerical code and the results 
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of few test problems. The results of the numerical simulations of the magnetized protostellar 
cloud collapse is presented in Section 6. In the conclusion (Section 7) we discuss briefly the 
main properties of the presented numerical code and general perspectives of its development 
in future. 



2 Problem statement 

2.1 Basic equations 



The system of the gravitational MHD reads: 

9P . ^ 



^ + V,{pv') = 0, (1) 
^(p^O + ^kipv'v' + g^'P - a'') = -pV^$ , (2) 
^-[V,[v,B]]^ = 0, (3) 



dt V^^ + + + + 7 + T ) + ^] 



^-fw'^Vk^, (4) 
V^$ = 47rG'p, (5) 



where 



is the Maxwellian tension tensor, g'^^ is the metric tensor. The other variables in this system 
(and further) have its usual physical sense. The system (1-6) is written in covariant form. 
Therefore it can be rewritten easy for the case of arbitrary curvilinear coordinates. 

To close the system (1-6) we use the equation of state of a perfect gas P = (7 — l)p£, where 
7 is the adiabatic index. The equation of energy (4) in the case of isothermal plasma should 
be excluded and the pressure should be determined from the relation P = c\^p, = RT/fj,, 
where ct is the isothermal sound speed, fi is the molecular weight. 

The collapse of a rotating magnetized protostellar cloud can be studied in the frame of 
axisymmetric approximation if the initial uniform magnetic field B is parallel to the vector 
of angular velocity fl. In this case we can use the cylindrical coordinates {r,(p,z), axis z 
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being parallel to the vectors B and fl. The variables in (1-6) do not depend on the azimuthal 
coordinate if, so equations (1-4) can be written then in divergent form as follows: 

du dF dG „ 

where the vectors of the conservative variables u, fluxes F, G and sources R are determined 
by expressions: 

u={rp, rpvr, r^pv^, rfw^, 

S,, S^, r5,, r (p£ + py + — j I , (8) 



F = 



An 



r pVrVz 



BrBz 

An 



, , VrB^ - v^Bz , r {VrBz - VzBr) , 



/ P v2\ B2 

rpVr \ e ^ H — + rvr- r— (v • B) 

\ p 1 I 47r 47r 



(9) 



G 



rpvz , r pVzVr - 



BzBr 

An 
B^ BV 



, r [pv^Vz 



B,pBz 
An 



r ypvl + B + — - ^ j , VzBr - VrBz , VzB^ - v^Bz , , 

/ P B^ E,, 

rpvz l^+- + yl+™^i;^-^i^(^-B) 



(10) 



?i-r— -r — 
Stt An ^^9r ' ' dz 



R= 0, pvl + P + — --^-rp—, 0, -rp— , 
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0, 0, 0, -rpVr-^-rpVz-^\ . (11) 
Here index T denotes the operation of the transposition. 

2.2 Initial and boundary conditions 



Let us consider the problem of initial and boundary condition on the example of the in- 
terstellar (protostellar) clouds collapse. Under the condition of axial and equatorial symme- 
tries we can considered the 2D computational domain in the cylindrical coordinate system: 
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D = {0<r<R,0<z<Z) with the characteristic sizes R and Z. We suppose that 
contracting cloud is located inside the computational domain and cloud boundary does not 
coincide with the grid boundary. 

The conditions on the external boundaries of the computational domain are defined by 
cloud contraction from the given volume. In this case the normal components of the velocity 
v„ on the external boundary are equal to zero. For density p, internal energy s, magnetic 
field B and tangential components of the velocity v,- we can write boundary conditions 
as d/dn = 0, where n is the outer normal unit vector. The internal boundary conditions 
correspond to the axial and equatorial symmetries {vr = — on the axe and d/dz — 0, 
Vz = = = on the equator). 

Initial conditions for the collapse of the rotating magnetic cloud are determined by values 
of the mass Mq, temperature Tq and parameters e^, e^,, e^^. These parameters are ratios of 
the internal, magnetic and kinetic energies of rotation to the absolute value of gravitational 
energy, respectively. The parameters e^, e^, allow us to obtain the following similarity 
relations for the variables characterizing the initial cloud state: 

i^ = 8.9in.l0^«.(^)(^)-'c.. (12) 
Po^ 6.7475.10-V(^)"(^)%-cm-, (13) 
..^ 7.5204. 10-V-.-(^)-(^)''%-.. (14) 

£.0 = 1.2342. 10-V4'^(|i)-'(^)=G. (15) 

These similarity relations allow to generalize the results of simulation with the specified 
values et, e^, on the contracting clouds with other values of mass and temperature. 

The minimal initial values of the parameters e^, e^, depend on the critical mass M*. We 
can obtain the value of critical mass from virial theorem as follows: 

3 / 3 

= + -(7 - + Y (^'^ + 2 ~ + ■ 

3 Finite-difference scheme 

3.1 TVD schemes 

Let us construct the finite-difference scheme for adiabatical and isothermal MHD fiows 
using the total variation diminishing (TVD) principle. TVD principle for schemes approxi- 



5 



mating the linear hyperbolic system is an extension of monotonia principle for the schemes 
approximating the linear advection equation. 

The general idea of this principle consists in the numerical implementation of the total 
variation diminishing for the finite-difference Riemann invariants [23]. An efficient way of 
the monotonic schemes construction was proposed by Vyaznikov et al. [24] . It consists in the 
following steps: 

(1) The selection of a monotonic 'base' scheme of the first order of approximation. 

(2) The transformation of the 'base' scheme to the high resolution scheme by addition of 
the antidiffusional terms. 

(3) The construction of the antidiffusional coefficients for monotonicity guarantee. 

The constructed TVD scheme should have a high resolution in the smooth regions of the 
solution but the order of approximation can decrease in the regions of large gradients. 



3.2 Advection equation 



To investigate the general ideology of the TVD scheme construction and its basic properties 
let us consider the linear advection equation: 

^ + ^^ll = , const , (16) 

that is defined in the domain D — {t > 0, —oo < x < +00} with the initial values p{x, 0) = 
Po{x). We introduce an uniform spatial grid {xi = ih, i = 0, ±1,±2, ...} with the fixed 
stepsize h and consider the following class of 'base' schemes for the equation (16): 



4 h h 

The parameter w is a coefficient of numerical viscosity, while the coefficients Q;j+i/2 > 
determine the antidiffusional terms. 

If all Q;i+i/2 = the scheme (17) gives the 'base' scheme of the first order of approximation 
in space and time. The condition of monotonicity of this scheme is w > [ul. In the case 
w — \v\ the 'base' scheme transforms into the classical donor cell scheme [25]. 

To obtain the conditions of monotonicity in the case of Q;i+i/2 7^ we have to use the 
principle of maximum [26], that for the scheme (17) gives the inequalities: 

1 _ , 1 «r-i/2 ^ ^ 

^ ^ -"-1-1/2 
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i<V. + 5il^S0, (18) 

^ ^ -"-1+1/2 

\ , / 1 1 \ 



where 



l«i-l/2\ w + v ( 1 + la»>i/2 
^ A T 1 (yT , ,„ : — '— 



1 + 1<V2^ 

<1, 



-"-1+1/2 — n _ „n ' -"-i-1/2 



Pi+l - Pi ' Pi - Pi-1 

are smoothness analyzers. Note that last inequality (18) determines also the stability condi- 
tion for the scheme (17). 

To monotonize the scheme (17) a number of approaches were suggested (see, e.g., [27-30]). 
Vyaznikov et al. [24] proposed to consider the coefficients a as a piecewise-linear functions of 
the analyzers R. These functions must satisfy to the following conditions: 1) a = 0, R < 0; 
2) < a < 2, R > 0. Let us take a{R) in the form of two-parametric set of piecewise-linear 
functions proposed by Chakravarthy and Osher [31]: 

<^{R) — —TT^ ■ P ■ minmod(l//9, R) + — r-^ • minmod(/9, R) , 



where 



minmod(x,y) = ^(sign(x) sign(y)) min(]x|, |y|) 



The function minmod(x, y) equals to zero if x and y have the different signs and to argument 
with the minimal modulus otherwise. The function a{R) is linear in the region of the smooth 
solution: l//3<R<p, 

a{R) = — ^ + —^R (19) 



and increases from value ai = ((1 + il))(3 + {I - ip))/{2p) to a2 = ((1 + ip) + {1 - ip)(3)/2. 

The scheme (17) has the second or third order of approximation on space in the regions 
of smooth solution. This result can be obtained by of expanding of the scheme (17) to the 
Taylor series in the neighborhood of point M — (t — t^,x — Xi). The corresponding first 
difference approximation for the scheme (17) is: 

Therefore at the value ip = 1/3 in the region of smooth solution the scheme (17) attains the 
third order of approximation in space. For the other values of -0 the scheme (17) has the 
second order of approximation. 
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The parameters ai and 0.2 determine the range of values of the function q:{R) in the 
smoothness region (where i? ~ 1). The optimal value of the parameter (3 attains at 02 = 2 
when /? = Prnax = (3 — V')/(l — V') • In this case we can find that ai = 2/(3 — -0), q;2 = 2. The 
stability condition is determined by last inequality (18). It can be rewritten in the following 
form: 



wr (w — v)r fa] (w — v)t . , , 

— I 7-; max < — > --; mm \ a Y 

h Ah R \R) Ah R ^ ^ 

iw-\-v)T f a 1 iw + v)t . . 
+ — 71-^ max <^ - ^ - ^ — min {«}<!. 
Ah R \R) Ah R ^ ^ - 

Prom this inequality and from the definition of a{R) we can obtain the stability condition 
of Courant, Priedrichs and Lewy (CFL condition, [32]): 

WT A 



In the case of optimal values of the parameters ■0 = 1/3 and (3 — A we can rewrite (20) by 
the following way: 

^<0.4. 
h 



The developed scheme can be rewritten more briefly in the conservative form: 

Pi^^ ~ Pi , /i+l/2 - fi-l/2 _ „ X 

T + h 



where 



1 — ip 

fi+i/2 = f-+i/2 ^ minmod(/^3/2, /^/i^i/s) 

- minmod(i^- 1/2> Pfz~+3/2) 

+ minmod(/+ 1/2, Pft-1/2) 

+ ^ minmod(/+ ,/2, Pfti/2) , (22) 

.0 _ fi+l + fi _ ^ fn _ n 

Ji+1/2 — 2 2 V ^* 

fi+1/2 ~ fi+1/2 ~ fi -I fi+1/2 ~ fi+'^ ~ fi+1/2 



and fluxes fi — vpf. 
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3.3 System of hyperbolic equations 



The developed methodology of difference scheme construction for the linear advection 
equation can be generalized onto the system of the linear hyperbolic equations: 

where u is the vector of the unknown variables and A is the matrix of the constant coefficients 
(hyperbolicity matrix). This system can be rewritten in the conservative form: 

du <9F _ ^ 

dt dx 

with the vector of fluxes: F = Au. 

The system of equations (23) can be represented in the equivalent form of the equations 
for the Riemann invariants 

S" = {r ■ u) (24) 

as 

where 1" are the left eigenvectors of the matrix A and Aq are its eigenvalues. The original 
variables u can be obtained using the inverse transformation: 

TV 

n^J2^aS\ (26) 

a=l 

where are the right eigenvectors of the matrix A. 
We approximate each equation in the system (25) by the difference scheme (21): 

Qa,n+1 Qa,n zra zra 

/ 1 



where fluxes 



- minmod (//";_- /2, 

+ minmod(i/f/,/„ +/,) 
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ffO' _|_ 77" 7,, 

TJCifi i+1 ' i Q / Qcv.n QQ,n\ rra \ Qa,n ^ I \ I 

tra— r7-o:,0 rra H"<^+ H"C« r7-Q;,0 

-"i+1/2 — -"i+1/2 -"i ' -"i+1/2 — -"i+1/2 ■ 

To construct the scheme for the variables u we have to perform the transformation (24) in 
(27). Finally the scheme reads: 

^'^i+ll^^lizlll ^ Q ^ (28) 
T h 

where fluxes are: 

= - ^ E minmod(F-3/2, (5^^^,,,) 

- E minmod(F^- /^Ff-g/J 

+ E minmod(F^+ /3Ft+ /J 

^ a 

+ ^ E minmod(Ft+ /2, /9F^+ / J , (29) 

F°+V2 = ^^^^ - ^RWL (ur^, - <) , 

W = diag(wi,W2, . . . ,WAr) , 
FS/2 = ^ ± ^") (l. ■ - <)) r" , 

where L, R - matrixes of the left and of the right eigenvectors of the matrix A respectively. 

The difference scheme (28,29) satisfies to the total variation diminishing (TVD) principle 
[23] , since it is constructed on the basis of the monotonic schemes for the Riemann invariants. 
In the present paper we consider the Lax-Friedrichs schemes with the coefficients Wa — w > 
max{|Aa|}. For this case the fluxes in the base scheme reads [33,34]: 

F°,V2 = - ? - <) • (30) 



During the calculation of fluxes (29) the numerical procedure of the minmod-evaluation 
should be called very often. To decrease the computation time we consider a more simple 
scheme in which the additional fluxes are summed up. Finally the scheme (28) is determined 
by the fluxes: 
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where 



I — jh 

Fi+i/2 = F°+i/2 — minmod(Fr^3/2, PF;^i/2) 

- mmmod(F-^i/2> 

+ l±t minmod(F+ P^ty^) 

+ ^ mmmod(F+ ,/2, PK1/2) > (31) 



Tji— rpO TP Tri+ TP TpO 

* i+1/2 - * i+1/2 - * i+1/2 - * i+1 - -I^ i+l/2 



and the fluxes F°_^^/2 determined by the expression (30). 
This methodology can be generahzed easy onto the system of the nonhnear equations: 

du dF(u) , , 

The finite-difi^erence scheme (28) for the system (32) is determined by the same fluxes 
(30,31), but the coefficients w we should take in the form: 

w = = 0max{|Aa,i|, |Ac,,j+i|} , 0>1, (33) 

where are the eigenvalues of the hyperbolicity matrix A = dF/dn. 



3.4 MHD modification of the finite- difference scheme 



We construct now the finite- difference TVD scheme for the isothermal or adiabatical ID 
MHD in Cartesian coordinates. Wc assume that the magnetic field and the velocity have all 
three components: B = {B^, By, B^}, v = {vx,Vy,Vz}- 

The MHD equations in the conservative form, 

du dF _^ 
dt dx 

are determined by the vector of the conservative variables: 



r v2 en 

u = p , pvx: pvy, pv^, B^, By, B^, ^^^^Y^Snl 



T 



and by the vector of the fluxes: 
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BxBz 



, v^By - VyB^ , v^B^ - v^B^ , 

( P 'B^\ B 

pv^ f £ + — + y + — j - ^ {v^B^ + VyBy + V^B^) 



For the adiabatic case the hyperbohcity matrix A is 
/ 







1 







Brr 



A2I (3 - -f)Vx (1 - -f)Vy (1 - 7)^;^ ^26 ^27 



-VxVz Vz 







^61 -TT 



V ^81 A 



-f 







B 



82 ^83 



-Vy Vx 







^ 



^84 ^85 ^86 -487 '^Vx j 



where: 



A 2 2 

^21 = - 



A2a=(2-7)§, 



^27 = (2 - 7) 



An ' 



161 



VzBx - VxBz 



171 



^81 = -Vx + " 2) ^ ^ "^^^ ■ ' 

A, = + — + - - (7 - 1)^2 ^ 



7-1 2 

^83 = (1 - 7)'t^a;^^y " O^xtty , AsA = (1 " ^)VxVz " a^O^ , 
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Wc use here the sound speed, c = \J'^P/p and a = B/v^47rp. The calculations give the 
following expressions for the eigenvalues of the hyperbolicity matrix (34): 



Ai = , A2 = f X , A3,4 = f ± a. 



A, 



'5,6,7,8 




c2 + a2 1 
2 ^2 




(35) 



The zero value of Ai is the consequence of zero value of the flux F component corresponding 
to Bx- The value of A2 corresponds to the entropic wave, A3 4 correspond to the Alfven 
waves, and A5_6,7,8 correspond to the fast and slow magnetosonic waves. In the isothermal 
case the entropic wave disappears and instead adiabatical sound speed it should be used the 
isothermal sound speed ct- 

We should construct now the calculation rule for the coefficients w in (33). The values 
of these coefficients should be not less than the maximal modulus of all eigenvalues of the 
hyperbolicity matrix. It is clear that this condition can be satisfied by the following way: 



Therefore we can keep in the relation (33) only eigenvalues corresponding to the fast MHD 
waves. This allow us to adapt the scheme (28,30,31,36) to the simulation of the ID MHD 
flows. If the magnetic field is absent the fast magnetosonic speed equals the sound speed 
c, therefore the transition to the gasdynamics occurs correctly. With the help of the splitting 
technique we can construct the scheme for multidimensional (2D and 3D) MHD flows. One 
of this scheme for the axisymmetric case we will describe in the next section. 

4 Numerical method 

4..1 General scheme 

The computational grid is determined by points (r^, Zj) indexed by integer indexes i — 
0, 1, . . . , /; j = 0, 1, . . . , J. These points denote the centers of the grid cells. The coordinates 



Wi+i/2 = 0max{|i;^,j| + , \v^,i+i\ + '^x,i+i} , > 1 , 



(36) 



where 
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of the cell boundaries are defined as the arithmetic mean of corresponding coordinates of 
the cell centers and have half-integer indexes. The computational variables are related to 
the centers of the cells and the numerical fluxes are related to the cell boundaries. Above we 
developed the TVD scheme in the frame of ID MHD approximation. Now we generalize this 
approach onto the 2D MHD approximation. The 'base' difference scheme for equations (7) 
is the following one: 



, Fi+l/2,j - Fj-l/2j Gi,j+l/2 - Gi,j-l/2 _ -p 

At ^ Ar ^ Az ~ ' 



,3 



The main problem here is the definition of the calculation rule for the numerical fluxes across 
the cell boundaries. 

Let us consider for simplicity only the radial direction because of the expressions of fluxes 
for the vertical direction can be derived by similar way. The radial fluxes Fi+i/2j and Fi_i/2j 
should be calculated with the help of the formulae (30,31,33). Relatively determination of 
the coefficients w we should note the following. 

In the definition of the vectors u, F and G in the form (8-10) the radial coordinate r is 
present. Therefore to calculate eigenvalues of the hyperbolicity matrices we can close formally 
the system (7) by the following obvious equation: 

It can be shown that this extension of the original MHD system not changes the standard 
set of the eigenvalues. 

In fact we consider the following system of MHD equations: 

dt dr dz 

with the vector w determined as 

w = < r , p, pvr, pv^, pv^, Br, , B^, ^ ^ T J ^ Stt I ■ 

Obviously in Cartesian coordinates the eigenvalues of matrices A and B are determined by 
the same expressions (35). Let us consider another vector of the variables u that consists of 
r and eight components of the vector (8). Making the transformation of the equation (37) 
to the new variables u we obtain: 

du ^/"^^ B'*^^ R' 
dt dr dz 

where 

A' = T-^AT, B' = T-^BT, R' = T-^R, (38) 
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and matrix T = du/dw is Jacobian of transformation. The determinant of this matrix 
det T = 1/r is not equal to zero. Hence the formulae (38) describes the similarity transfor- 
mations of matrices A and B that does not change their eigenvalues. 

This discussion allows us to take the values of the coefficients Wj+1/2 in (33) in the following 
form: 



w. 



1+1/2 = (l)max{\Vr,i\ + *r,j , \Vr,i+l \ + *r,i+l}, > 1 , (39) 



where 



c2 + a2 + ,/(c2 + a2)2 - Aay 
^ ^ — ■ (40) 



These coefficients are determined by eigenvalues that correspond to the fast magnetosonic 
wave. The expressions (39) can be apphed both to the adiabatical and to the isothermal 
plasma flows. The energy equation should be excluded in last case and as c in (40) should 
be used the isothermal sound speed ct- 

Since the fluxes related to the different space directions are constructed independently, the 
^-part of the numerical scheme can be developed by similar way. 



4-2 Poisson solver 



Poisson equation for the gravitational potential is solved by ADI-method (alternating di- 
rections implicit method) of Douglas and Rachford [35-37] : 



eKr^ + (1 - e)A,,^p - ^ 



= (1 - e)Arr^ + eK,^^+' - P- 



where 9 is the parameter of scheme implicity, A^^j ^zz are second order finite-difference 
approximations of the corresponding parts of Laplacian operator in r- and ^-directions, 
respectively, Tp is the iteration parameter. 

Boundary conditions for the gravitational potential on the inner boundaries are regular: 
9$/9n = 0. Boundary condition on the external boundary is = ^ext where <^int is 
the computing potential, (^ext is the external potential of the protostellar cloud. It can be 
obtained from Poisson integral transformation of equation (5) solution: 
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where X = zij-"^ + z^) and Pi{x) is the Legendre polynomial of lih. order. All terms with 
odd / are absent in this representation as a result of the axial symmetry. 

4-3 Implementation of the divB = condition 

Another difficulty of the numerical simulation of the multidimensional MHD flows is con- 
nected with the necessity of the numerical implementation of the vanishing divergence con- 
dition for the magnetic field 

divB = 0. (41) 

Numerical errors due to discretization and finite-difference differential operator splitting 
can lead to nonzero divergence over the time. This problem can be solved by many ways [38]. 
In the first approach the 8-wave formulation of the MHD equations is used (see also [11]). 
Another way uses the constrained transport suggested by Evans and Hawley [39]. Finally, 
the third approach is the projection scheme [40]. Our numerical scheme uses the last method 
(see [18,41]). 

Denote as B* the magnetic field obtained by the numerical solution of the induction equa- 
tion (3). This 'numerical' magnetic field contains the numerical error b, that consist of vortex 
and potential parts: 

b = [V, a] + V0 . 

Therefore exact solution of induction equation is: 

B = B* - b . (42) 

Vortex part of numerical error [V, a] does not contributes to divergence of the magnetic field 
B. Substitution of (42) into (41) shows that potential (j) must satisfy the Poisson equation 

V^0 = divB*. (43) 

Resolving this equation numerically for the defined 'numerical' magnetic field B* for each 
time step we calculate the divergence free magnetic field with the help of formula 

B" = B* - V0. 

Further we will assume that B = B''. The algorithm is implemented in our numerical code 
as follows. Denote the central cell with indexes (i, j) by the index C for brevity (see. Fig. 1), 
the neighbor cells by indexes W, N, E, S and intermediate coordinates (cell boundaries) by 
double indexes CE, CS etc. 
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Fig. 1. A typical grid cell 
Then we can approximate the divergence of magnetic field by the expression: 
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(44) 



where intermediate values are calculated by the formulae 



etc. 



The finite-difference approximation of the Laplacian operator is 
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It is easy to see that if the potential (f) is the exact solution of the finite-difference Poisson 
equation (43) then the magnetic field B (for which the finite-difference divergence (44) is 
equal to zero exactly) can be obtained by the formulae: 
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The Poisson equation (43) for the potential (f) is solved in our scheme by the ADI-method 
that is described in the previous section. As the additional boundary conditions on the inner 
boundaries (r = 0, 2; = 0) we use the relations dcp/du = 0. On the outer boundaries we use 
the boundary condition = since the external field has identically vanishing divergence. 



5 Test computations 

5. 1 Background 

The numerical code 'Moon' was developed on the basis of the proposed finite- difference 
scheme. It can simulate the ID and 2D self gravitating MHD fiows. The code is implemented 
by the programming language C-|— I- and completely object-oriented. Source text of the pro- 
gram uses only standard C-|— I- elements and does not relate to concrete compiler. The kernel 
of the code uses intensively the object-oriented library 'Numerical Tool Box' (NTB 3.5). 
This library was developed by A.G.Zhilkin in 1996 for more easy creation of the complex 
numerical codes. To check the properties of our code we tested it on the number of test 
problems with known (exact or approximate) analytical solutions. 

5.2 One- dimensional tests 

5.2.1 Test 1. Linear advection 

The linear advection of the square density profile was simulated to check the phase-error 
properties of the numerical scheme. Figure 2 gives the density profile after 100th time step 
without (triangles) and with fiux correction (dots). The computational grid consists of 100 
cells. The analytical solution is drawn by the sohd line. 

It is seen that the Ist-ordcr of approximation method has a large numerical diffusion. 
This diffusion smears the initial profile rapidly, the width of smearing being increased during 
calculation. The method of the 3rd-order of approximation produces a profile that is smeared 
onto 3-4 cells only and conserves its form during a long time. We can say therefore that 
amplitude and phase errors in the scheme of the 3rd-order of approximation are selfconsistent 
and its order of approximation is not less actually then the approximation order of the so- 
called LPE-schemes (little phase error) (see [6]). In our opinion this circumstance is connected 
with the successful selection of the antidiffusion limiters (22) in the scheme (21). 

5.2.2 Test 2. Decay of arbitrary gasdynamical discontinuity 

Numerical simulation of the arbitrary discontinuity decay (Riemann problem) allows us 
to check the scheme properties relating to the resolution of the shock waves, contact dis- 
continuities and rarefaction waves. At the initial time moment the computational domain 
along the x-axis was divided on two subdomains A and B. In the initial state the pressure 
and density have the following values: in subdomain A — = 1, = 1; in subdomain B 
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Fig. 2. Linear advection of squared density profile. The density distributions is shown after 100 
time-step computed by the scheme of the Ist-order of approximation (triangles) and by the scheme 
of the 3rd-order of approximation (dots). The solid line corresponds to the analytical solution 

— Pb = 0.125, Pb = 0.1. The velocity in both domains va = vb = ^ and adiabatic index 
7 = 5/3. Figure 3 shows the analytical solution (sohd line) and the results of the numerical 
computation (dots) at the time moment t — 0.164. As a result of the discontinuity decay the 
rarefaction wave 1 propagates to the left, contact discontinuity 2 and shock wave 3 prop- 
agates to the right. The Figure shows that the scheme resolves the shock wave sufficiently 
well, but gives the small non-physical oscillations on the contact discontinuity. To smooth it 
is necessary to use additional hmiters of antidiffusion fluxes. 



5.2.3 Tests 3, 4- Decay of arbitrary MHD discontinuity 

The numerical simulation of an arbitrary MHD discontinuity decay is the unique possibility 
to check the MHD properties of the difference scheme. The problem statement in this test 
is similar to the problem statement in test 2. The initial parameters for this test and its 
exact analytical solution were taken from paper of Ryu and Jones [18] (variant la): 7 = 5/3, 

Pa = 1, Pb = 1, Pa = 20, Pb = 1, v^^a = 10, v^^^b = -10, By^A = 5, By^B = 5, = 5. 
The values of another variables are taken to be zero. This situation can arise, for example, 
at the collision of two gaseous masses moving forward each to other at the angle 7r/4 to the 
direction of the magnetic field. 

The results of the numerical computation and the analytical solution of this problem at 
the time t — 0.08 are shown on the Fig. 4. As a result of the discontinuity decay two fast 
MHD shock waves (1 and 5) arise and propagate in the opposite directions. The rotational 
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Fig. 3. Decay of an gasdynamical discontinuity. The density distribution (dots) is shown at the time 
moment t = 0.164. The sohd hne corresponds to the analytical solution. The digits correspond to: 
1 — rarefaction wave, 2 — contact discontinuity, 3 — shock wave. In the regions of smooth flow the 
dots is displayed more seldom (each 7th point). In the regions of contact discontinuity and shock 
wave all points of numerical solution are displayed 

(alfvenic) discontinuity 2, the contact discontinuity 3 and the slow MHD shock wave 4 are 
formed between this waves. The slow MHD shock wave propagates immediately after the 
fast MHD shock wave to the region of the gas with the larger pressure. 

In another variant of this test we consider the Riemann problem with the initial state vector 
{p,v^,Vy,v„By,B-„P) (see [18], variant 4c): region A — (0.65,0.667,-0.257,0,0.55,0,0.5), 
region 5 — (1, 0.4, -0.94, 0, 0, 0, 0.75) with B,, = 0.75 and 7 = 5/3. The switch-of slow shock 
wave is formed in this case. The results of computations at the time t = 0.15 are shown on 
Fig. 5. A fast MHD shock wave 1, switch-of slow MHD shock wave 2, contact discontinuity 
3 and gasdynamical shock wave 4 are formed. 

These tests show that the scheme approximation of fast and slow shock waves as well as 
the rotational (alfvenic) discontinuity is good enough. On the contact discontinuities and on 
the switch-off shock wave small oscillations occur. 

5.2.4 Test 5. Alfven wave 

The problem of propagation of the finite amplitude Alfven wave (the simple alfvenic wave) 
is interesting because it has the analytical solution (see [42]). The comparison of the numerical 
and the analytical velocity profiles and the magnetic field allows us to make conclusions about 
behavior of the phase and the amplitude errors in the MHD case. We solve the following 
problem. At the initial time moment the wave is concentrated in the region xl < x < xr, 
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Fig. 4. Decay of an MHD discontinuity (variant 1). The solution is corresponded to time moment 
t = 0.08. The digits correspond to: 1,5 — fast MHD shock waves, 2 — rotational (alfvenic) discon- 
tinuity, 3 — contact discontinuity, 4 — slow MHD shock wave 



where 



7-, ^ { xr — x 7r\ „ { xr — x 7r\ 

By = 10 sm TT , = 10 cos TT , 

\ xr-xl 2J \ xr-xl 2J 



By Bz 

Vy - — - - 



The wave propagates in the gas with parameters p = 1, P = 1, = Bx = 1^ = lA with 
velocity a = B^/ \/47rp without profile distortion. In wave region the tangential component 
of the magnetic field rotates with the conservation of its absolute value. Figure 6. shows the 
numerical and analytical distributions By and Bz at the time moment i = 0.3. 

The comparison of these curves with each other allows us to make the conclusion that even 
in the case of full MHD system the essential smearing of the velocity and the magnetic field 
profiles is absent during a sufficiently long time of the computation. Hence the phase errors 
balance the amplitude ones with the high order of accuracy. 



5.3 Two-dimensional tests 
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Fig. 5. Decay of an MHD discontinuity (variant 2). The solution is corresponded to time moment 
t = 0.15. The digits correspond to: 1 — fast MHD shock wave, 2 — slow MHD shock wave, 3 — 
contact discontinuity, 4 — gasdynamical shock wave 

5.3.1 Test 6. Two-dimensional advection 

To check the diffusion properties of the developed 2D code we obtain the numerical solution 
of the 2D advection equation (see similar computations of Munz [43]) 

in box domain (0 < a;, |/ < 1). The analytical solution of this equation for the velocity profile 

Vx ^ -{y - Voji^ , Vy^{x-XQ)UJ. 



describes the rotation of the density profile around the point (xq, Uq) with the angular velocity 
u). We use in this computation the following values: Xq = Hq = 0.5, a; = 1. Hence to the 
time t — 2'K the original density profile should carry out one full rotation and return to 
its initial position. The computation were done on the uniform grid with the cells number 
= 100 X 100. We compare the initial state with the profile that obtained after one full 
rotation. 

Figure 7a shows the original density profile. It is split on two regions. In the first region 
the density is equal to zero and in the second region the density is equal to 1. The boundary 
of the regions is circle with the cut-out rectangle. The center of the circle is not coincides 
with the center of the domain. Figure 7b shows the density profile at the time t — 27: after 
one full revolution. The profile (as in the ID case) is smeared onto 3-4 cells. Therefore we 
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Fig. 6. Alfven wave. The numerical (circles) and the analytical (solid lines) distributions of the 
transverse components of the magnetic field are shown at the time moment t = 0.3 




Fig. 7. Two-dimensional advection. The density distributions are shown at the initial time moment 
(a) and after one full revolution t = 2tt (b) 



can conclude that 2D variant of our scheme has the good selfconsistent amplitude-phase 
properties as well. 
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5.3.2 Test 7. Expansion of blast wave 

The test problem of blast wave expansion in the uniform medium is very crucial for 2D 

numerical code due to spherical symmetry of the wave. Simulation the blast wave expan- 
sion allows us to check how precisely our 2D axisymmetrical numerical code can resolve 
spherically-symmetrical problems. We compute, in particular, the expansion of the supernova 
remnant on the adiabatic stage. This stage is described by Sedov-Taylor solution [44-46], 
the approximation being true while the radiation cooling is small. 

We use the following initial parameters. The supernovae explosion energy is equal to 4 • 
10™ erg, the density of the interstellar medium po ~ 2 - 10~^^ g-cm~^, the spatial scale of the 
computational domain Rq k, 4.03- 10^^ cm, the pressure in the central region differs from the 
pressure in the external medium at the initial time moment approximately by eight order of 
magnitude. The scales of the computational values are equal to: for distance — 13 pc, for 
velocity — 8.9 x 10^ cm • s"^, for pressure — 1.1 x 10~^ erg • cm~^, for time — 1510 Yrs. The 
expansion of the supernova remnant can be followed numerically over the time t — 10600 Yrs. 
The radius of the shock wave increases from 3.6 pc to 12 pc during this time. 

Figure 8 shows the dependencies of the radius R, velocity U and the pressure P of the 
shock wave on the time. All curves correspond to power law: {R,U,P) oc t^. We obtained 
the following exponents: for radius of the shock wave kR = 0.403, for velocity ku = —0.589, 
for pressure kp — —1.139. The analytical values are kji — 0.4, ku = —0.6 and kp — —1.2. 
We see that the numerical values of the exponents are very close to exact Sedov-Taylor 's 
values. The condition of the spherical symmetry over the computations satisfies within the 
accuracy of 0.3%. 



5.3.3 Tests. Spherically-symmetrical free-fall collapse 

For simulation of free-fall (pressure-free) collapse as an initial state we take the uniform 
spherically-symmetrical cloud with the mass M = IMq and with the initial density po = 
10^1^ g ■ cm~^. Figure 9 shows the evolution of the numerical density (dots) of the cloud 
over the time in comparison with the analytical solution (solid line). We can see that the 
numerical points of density coincide with the high accuracy with values of the analytical 
solution. The cloud contracts to the point at the free-fall time tff — Y^37r/(32G'po)- The 
computation is carried out to time moment when the cloud size decreases to the size of one 
grid cell. The difference of numerical and analytical values of density is not greater than 1% 
at this time and the law of the spherical symmetry conservation satisfies within the accuracy 
of 0.4%. 



5.3.4 Test 9. The contraction of polytropic cloud with 7 = 4/3 

The simulation of the selfgravitating polytropic cloud dynamics with the adiabatic index 
7 = 4/3 allows us to check the quality and the accuracy of the code in the simulation of the 
selfgravitating MHD flows. The hydrostatic equilibrium of polytropic cloud with 7 = 4/3 
is indifferent one since the total energy, its flrst and second variations are equal to zero 
(see [47]). The density distribution of this cloud can be presented by analytical form as 
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Fig. 8. Expansion of blast wave. The solid lines correspond to the numerical dependencies of the 
radius R (curve 1), velocity U (curve 2) and pressure P (curve 3) of shock wave behind its front of 
the time. The corresponding analytical dependencies are shown by dotted lines 

p(r) = pcX^{r/R), where pc is the central density, R is the cloud radius, x{^) = (^si^/^i) 
{Osi^) is the Emden function of index 3 and ~ 6.9 is its first positive zero). 

To check the accuracy of the conservation of mass, momentum, and energy by our scheme we 
solve numerically the equations of the gravitational gasdynamics with the initial conditions 
corresponding to the specified hydrostatic state of the cloud. We computed several models 
with different numbers of cells of the spatial grid. Figure 10 shows the density profiles in 
the case of = 100 x 100. From the comparison of numerical (dots) and analytical (solid 
line) density profiles on this Figure we conclude that the relative error of the numerical 
solution compared with the analytical one at the time t = l-5t// is less than 2%. The law of 
total mass conservation over the computations satisfies within the accuracy of 0.003% and 
the law of the total energy conservation satisfies within the accuracy of 2%. 

This problem was solved on the grids with the number of cells = 20 x 20, 80 x 80 and 
120 X 120 as well. We can say that the relative error of the numerical solution obtained by the 
time t = l.St/j decreases with the increasing of the cells number. In the case of = 20 x 20 
it equals to 20% while in the case of A^ = 80 x 80 — 7%, and for A^ = 120 x 120 it is less 
than one percent. We can conclude therefore that the convergence of the numerical solution 
to the analytical one with the rise of the grid resolution takes place. 

In another variant of this test we specified the deviation from the hydrostatic equilibrium 
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Fig. 9. The free-fall collapse. The numerical (dots) and the analytical (solid line) dependencies of 
the cloud density of the time are shown. As a unit of density it is used its initial value. As a unit 
of time it is used the cloud free-fall time tff 
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Fig. 10. Hydrostatic evolution of the selfgravitating polytropic cloud with 7 = 4/3. The analytical 
(solid line) and the numerical (dots) profiles of the density are shown at the time moment t = 1.5t// 
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Fig. 11. Contraction of the polytropic selfgravitating cloud with 7 = 4/3. The distributions of 
the density (a) and the velocity (b) are shown for the different moments of time. The circles and 
triangles correspond to the numerical solution and the solid lines correspond to the analytical one. 
The digits 1,2,3 correspond to the moments of time O.Oto, 0.69to and 0.82to, respectively. 

at the initial time moment by the velocity distribution v{r,t) = —Ur/R. The cloud radius 
in this case evolves in time in accordance with the law R{t) — Rq — U -t and the evolution 
of the density and the velocity in the cloud can be described by the following analytical 
dependencies: 




To the time to — Rq/U the cloud will contract to the point. 

Figure 11 shows the numerical distributions of the density (circles) and the velocity (trian- 
gles) for different times. The solid lines correspond to analytical curves. Prom analysis of this 
Figure we can conclude that the numerical solution with the high accuracy coincides with the 
analytical one while the grid resolution is sufficient for the adequate spatial approximation 
of the central region. For the variant with = 100 x 100 this corresponds approximately to 
the time moment t — O.Sto (see Fig. 11). 



6 Collapse of magnetized protostellar clouds 

6. 1 Introduction 

Protostellar clouds as the cores of interstellar molecular clouds consist of neutral (atoms, 
molecules and cosmic dust particles) and charged (electrons, ions and charged dust) compo- 
nents. The equations of similar multicomponent mixture are very complex for the numerical 
simulations. 
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The more rough ideal MHD approximation can be used as a first approach. Contracting 
cloud is transparent to infrared self-radiation on the early stages of collapse. Therefore the 
initial phases of protostellar collapse can be considered as the isothermal ones with the good 
approximation. This approximation works while the optical depth of the collapsing cloud 
core becomes comparable with the unity. 

We performed the computation of the isothermal magnetized cloud collapse using the 
uniform grid with the cells number N — 120 x 120 and Courant number C = 0.3. We use 
the following values of parameters. The initial ratio of the thermal and kinetic rotational 
energies to the absolute value of the gravitational energy = 0.386, e^j = 0.0, respectively. 
The mass of the protostellar cloud Mq = l.SM© and its temperature Tq = 10 K. One the 
basis of the relations (12-15) we can conclude that these parameters correspond to the cloud 
with the initial density po — 5.18 • 10~^^ g • cm~^ and radius Rq — 5.17 • 10^^ cm. As a result 
of computations we calculate the flatness degree (ratio of the cloud thikness to its radius) 
and the central density up to the free-fall time. 

6.2 Kinematics 

The simplest case of the collapse corresponds to infinitesimally weak magnetic field. If the 
initial value < 0.5 ■ 10^"^ then we can neglect the influence of the magnetic fleld on the 
collapse dynamics. The magnetic fleld can be considered as a passive admixture and for the 
given velocity fleld and may be determined from the equation of induction (3). 

The results of computations show that in this kinematic approach the magnetic fleld B (x 
with k — 2/3 since of the cloud contraction is practically spherically-symmetrical one. Note 
that in the cloud envelope this law fulflUs approximately and in the forming core it satisfles 

exactly [48]. The computations conflrm the early obtained conclusion that the magnetic 
fleld in the cloud acquires with the time the quasi-radial geometry (see Fig. 12) [48]. This 
geometry of the fleld becomes distinct already to the free-fall time. 

6.3 Dynamics 

In the case of intermediate magnetic flelds (0.1 < < 0.5) the collapse picture drastically 
differs from the kinematic statement. The gas particles moving across the magnetic fleld 
lines are deviated to the equatorial plane by the electromagnetic forces. As a result the flat 
(disk-like) structure is formed. Figure 13 shows the typical density proflles and geometry of 
the magnetic fleld lines at the moment of time equal to the magnetic free-fall time 



(see [48]). The magnetic fleld in the disk is almost quasi- uniform and the ratio Br/Bz is 
small. 
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Fig. 12. The density distribution and the picture of the magnetic field fines in tfie central region at 
tfie time moment t = tff in tfie case of the kinematic cloud model 

Numerical simulations show that at the free-fall time the flatness degree of the disk is 
proportional to e^^"^ and the central density is inversely proportional to em- 



6.4 Quasi-steady-state contraction 



The growth of the initial magnetic field leads to the increase of the electromagnetic forces 
that support the cloud from the contraction. Therefore at the large initial values e^, the cloud 
relaxes to the quasi-steady-state equihbrium. Such a cloud should evolve within the diffusion 
time scale that is determined by the ambipolar diffusion of the magnetic field. The cloud 
loses the magnetic flux by the ambipolar diffusion and slowly contracts retaining the quasi- 
hydrostatic equilibrium. As a result this evolution the non-uniform profile of the density is 
formed. The lifetime of such molecular cloud cores can reach 10^ — 10^ years. 



7 Conclusion 



The explicit conservative finite-difference TVD scheme of high resolution for the solution 
of various MHD problems is constructed. The numerical MHD code 'Moon' on the basis of 
this scheme is developed that can simulate the MHD flows in the ID and 2D approximations. 
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Fig. 13. The density distribution and lines of equal magnetic flux in the central region of the cloud 
with em = 0.1 at the time moment t = tfm 

The results of extensive tests show that the code is well adapted to the simulation of many 
MHD problems of plasma physics and astrophysics. 

To check the dissipative properties of the developed scheme the problems of ID and 2D 
advection are solved. The results of these tests show that the scheme smears the initial 
density profile onto 3-4 cells. In the test computation of the blast wave the condition of the 
spherical symmetry satisfies within the accuracy of 0.3%. 

To analyze the convergence of the numerical solution to the analytical one the problem of 
the hydrostatic equilibrium of the selfgravitating polytropic cloud with the adiabatic index 
7 = 4/3 is computed on the grids with various number of cells. The obtained results allow 
us to conclude that the errors of the numerical solution decrease to the permissible values on 
the grids with the number of cells not less than N = 100 x 100. In the case > 120 x 120 
the relative error of the numerical solution at the time moment t = 1.55t// is less that 1%. 
The numerical solution in this case coincides with the good accuracy with the analytical one 
while the grid resolution in the central region is sufficient for acceptable approximation of 
the grid values. 

The numerical computations of the collapse of the magnetized protostellar cloud confirms 
the earlier obtained results in the framework of 1.5D approximation (see [48]). The magnetic 
field in the cloud with the initially small intensity acquires the quasi-radial geometry after 
some time. The intermediate magnetic field leads to flattening of the collapsing clouds on 
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the late stages of contraction. For the strong initial magnetic field the cloud relaxes to quasi- 
hydrostatic equilibrium. Such a cloud must evolve within the difi^usion time scale. 

Note that we can not simulate numerically the advanced stages of the collapse because the 
described code 'Moon' uses the uniform grid only. But investigation of these stages is very 
important and interesting problem to understand the physics of the formation and further 
evolution of young stellar objects. We are elaborating now the adaptive mesh refinement 
(AMR) numerical algorithm for the simulation of advanced stages of rotating magnetized 
protostellar cloud collapse. The AMR-approach uses the tree-threaded collection of subgrids 
with the subsequently increasing grid resolution. We are constructing a new 2D MHD numer- 
ical code 'Megalion' that based on the described in this paper high-resolution TVD scheme 
for the gravitational compressible MHD fiows. This code will allow us to simulate various 
selfgravitating MHD fiows with the large gradients and fast time variations. New numer- 
ical code will take into account also the processes of ambipolar diffusion, non-stationary 
ionization and heating/cooling. 
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